COVID-19 Analytical Report
1. Data Preprocessing
1.1 Gloabl Confirmed Cases and Deaths Time-Series by Country
ts.confirmed = covid19.data("ts-confirmed")## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## -------------------------------------------------------------------------------- ts.confirmed = ts.confirmed %>%
group_by(Country.Region) %>%
select(-c(Province.State, Lat, Long)) %>%
summarise_all(funs(sum)) %>%
arrange(desc(.[[ncol(ts.confirmed)-3]]))## Warning: funs() is soft deprecated as of dplyr 0.8.0
## Please use a list of either functions or lambdas:
##
## # Simple named list:
## list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`:
## tibble::lst(mean, median)
##
## # Using lambdas
## list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once per session. head(ts.confirmed)## # A tibble: 6 x 97
## Country.Region `2020-01-22` `2020-01-23` `2020-01-24` `2020-01-25`
## <fct> <int> <int> <int> <int>
## 1 US 1 1 2 2
## 2 Spain 0 0 0 0
## 3 Italy 0 0 0 0
## 4 France 0 0 2 3
## 5 Germany 0 0 0 0
## 6 United Kingdom 0 0 0 0
## # … with 92 more variables: `2020-01-26` <int>, `2020-01-27` <int>,
## # `2020-01-28` <int>, `2020-01-29` <int>, `2020-01-30` <int>,
## # `2020-01-31` <int>, `2020-02-01` <int>, `2020-02-02` <int>,
## # `2020-02-03` <int>, `2020-02-04` <int>, `2020-02-05` <int>,
## # `2020-02-06` <int>, `2020-02-07` <int>, `2020-02-08` <int>,
## # `2020-02-09` <int>, `2020-02-10` <int>, `2020-02-11` <int>,
## # `2020-02-12` <int>, `2020-02-13` <int>, `2020-02-14` <int>,
## # `2020-02-15` <int>, `2020-02-16` <int>, `2020-02-17` <int>,
## # `2020-02-18` <int>, `2020-02-19` <int>, `2020-02-20` <int>,
## # `2020-02-21` <int>, `2020-02-22` <int>, `2020-02-23` <int>,
## # `2020-02-24` <int>, `2020-02-25` <int>, `2020-02-26` <int>,
## # `2020-02-27` <int>, `2020-02-28` <int>, `2020-02-29` <int>,
## # `2020-03-01` <int>, `2020-03-02` <int>, `2020-03-03` <int>,
## # `2020-03-04` <int>, `2020-03-05` <int>, `2020-03-06` <int>,
## # `2020-03-07` <int>, `2020-03-08` <int>, `2020-03-09` <int>,
## # `2020-03-10` <int>, `2020-03-11` <int>, `2020-03-12` <int>,
## # `2020-03-13` <int>, `2020-03-14` <int>, `2020-03-15` <int>,
## # `2020-03-16` <int>, `2020-03-17` <int>, `2020-03-18` <int>,
## # `2020-03-19` <int>, `2020-03-20` <int>, `2020-03-21` <int>,
## # `2020-03-22` <int>, `2020-03-23` <int>, `2020-03-24` <int>,
## # `2020-03-25` <int>, `2020-03-26` <int>, `2020-03-27` <int>,
## # `2020-03-28` <int>, `2020-03-29` <int>, `2020-03-30` <int>,
## # `2020-03-31` <int>, `2020-04-01` <int>, `2020-04-02` <int>,
## # `2020-04-03` <int>, `2020-04-04` <int>, `2020-04-05` <int>,
## # `2020-04-06` <int>, `2020-04-07` <int>, `2020-04-08` <int>,
## # `2020-04-09` <int>, `2020-04-10` <int>, `2020-04-11` <int>,
## # `2020-04-12` <int>, `2020-04-13` <int>, `2020-04-14` <int>,
## # `2020-04-15` <int>, `2020-04-16` <int>, `2020-04-17` <int>,
## # `2020-04-18` <int>, `2020-04-19` <int>, `2020-04-20` <int>,
## # `2020-04-21` <int>, `2020-04-22` <int>, `2020-04-23` <int>,
## # `2020-04-24` <int>, `2020-04-25` <int>, `2020-04-26` <int>ts.deaths = covid19.data("ts-deaths")## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## --------------------------------------------------------------------------------ts.deaths = ts.deaths %>%
group_by(Country.Region) %>%
select(-c(Province.State, Lat, Long)) %>%
summarise_all(funs(sum)) %>%
arrange(desc(.[[ncol(ts.deaths)-3]]))
head(ts.deaths)## # A tibble: 6 x 97
## Country.Region `2020-01-22` `2020-01-23` `2020-01-24` `2020-01-25`
## <fct> <int> <int> <int> <int>
## 1 US 0 0 0 0
## 2 Italy 0 0 0 0
## 3 Spain 0 0 0 0
## 4 France 0 0 0 0
## 5 United Kingdom 0 0 0 0
## 6 Belgium 0 0 0 0
## # … with 92 more variables: `2020-01-26` <int>, `2020-01-27` <int>,
## # `2020-01-28` <int>, `2020-01-29` <int>, `2020-01-30` <int>,
## # `2020-01-31` <int>, `2020-02-01` <int>, `2020-02-02` <int>,
## # `2020-02-03` <int>, `2020-02-04` <int>, `2020-02-05` <int>,
## # `2020-02-06` <int>, `2020-02-07` <int>, `2020-02-08` <int>,
## # `2020-02-09` <int>, `2020-02-10` <int>, `2020-02-11` <int>,
## # `2020-02-12` <int>, `2020-02-13` <int>, `2020-02-14` <int>,
## # `2020-02-15` <int>, `2020-02-16` <int>, `2020-02-17` <int>,
## # `2020-02-18` <int>, `2020-02-19` <int>, `2020-02-20` <int>,
## # `2020-02-21` <int>, `2020-02-22` <int>, `2020-02-23` <int>,
## # `2020-02-24` <int>, `2020-02-25` <int>, `2020-02-26` <int>,
## # `2020-02-27` <int>, `2020-02-28` <int>, `2020-02-29` <int>,
## # `2020-03-01` <int>, `2020-03-02` <int>, `2020-03-03` <int>,
## # `2020-03-04` <int>, `2020-03-05` <int>, `2020-03-06` <int>,
## # `2020-03-07` <int>, `2020-03-08` <int>, `2020-03-09` <int>,
## # `2020-03-10` <int>, `2020-03-11` <int>, `2020-03-12` <int>,
## # `2020-03-13` <int>, `2020-03-14` <int>, `2020-03-15` <int>,
## # `2020-03-16` <int>, `2020-03-17` <int>, `2020-03-18` <int>,
## # `2020-03-19` <int>, `2020-03-20` <int>, `2020-03-21` <int>,
## # `2020-03-22` <int>, `2020-03-23` <int>, `2020-03-24` <int>,
## # `2020-03-25` <int>, `2020-03-26` <int>, `2020-03-27` <int>,
## # `2020-03-28` <int>, `2020-03-29` <int>, `2020-03-30` <int>,
## # `2020-03-31` <int>, `2020-04-01` <int>, `2020-04-02` <int>,
## # `2020-04-03` <int>, `2020-04-04` <int>, `2020-04-05` <int>,
## # `2020-04-06` <int>, `2020-04-07` <int>, `2020-04-08` <int>,
## # `2020-04-09` <int>, `2020-04-10` <int>, `2020-04-11` <int>,
## # `2020-04-12` <int>, `2020-04-13` <int>, `2020-04-14` <int>,
## # `2020-04-15` <int>, `2020-04-16` <int>, `2020-04-17` <int>,
## # `2020-04-18` <int>, `2020-04-19` <int>, `2020-04-20` <int>,
## # `2020-04-21` <int>, `2020-04-22` <int>, `2020-04-23` <int>,
## # `2020-04-24` <int>, `2020-04-25` <int>, `2020-04-26` <int>1.2 U.S. Confirmed Cases and Deaths Time-Series By State
nCov = load_nCov2019()
df.nCov = nCov$province
ts.us = df.nCov %>% filter(country == "United States") %>%
select("Date" = time, "State" = province, "Confirmed" = cum_confirm, "Deaths" = cum_dead)
ts.us$State = sapply(ts.us$State, Caps)
head(ts.us)## Date State Confirmed Deaths
## 1 2020-03-15 Washington State 769 42
## 2 2020-03-15 New York State 746 6
## 3 2020-03-15 California 431 6
## 4 2020-03-15 Massachusetts 164 0
## 5 2020-03-15 Florida 136 5
## 6 2020-03-15 Colorado 131 11.3 Global Aggregated Cases and Deaths by Country
agg.gloabl = covid19.data("aggregated")## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~agg.gloabl = agg.gloabl %>%
select(Country_Region, Confirmed, Deaths) %>%
group_by(Country_Region) %>%
summarise_all(funs(sum)) %>%
arrange(desc(Confirmed))
head(agg.gloabl)## # A tibble: 6 x 3
## Country_Region Confirmed Deaths
## <fct> <int> <int>
## 1 US 938154 53755
## 2 Spain 223759 22902
## 3 Italy 195351 26384
## 4 France 161644 22648
## 5 Germany 156513 5877
## 6 United Kingdom 149569 203811.4 U.S. Aggregated Cases and Deaths By State
agg.us = covid19.data("aggregated")## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~agg.us = agg.us %>%
filter(Country_Region == "US" & Province_State != "Recovered") %>%
select(Province_State, Confirmed, Deaths) %>%
group_by(Province_State) %>%
summarise_all(funs(sum)) %>%
arrange(desc(Confirmed))
head(agg.us)## # A tibble: 6 x 3
## Province_State Confirmed Deaths
## <fct> <int> <int>
## 1 New York 282143 22009
## 2 New Jersey 105498 5914
## 3 Massachusetts 53348 2730
## 4 California 42368 1689
## 5 Illinois 41777 1875
## 6 Pennsylvania 41153 17931.5 U.S. Active Cases, Deaths and Recovered Time-Series
ts.us.all = covid19.data("ts-ALL") %>%
filter(Country.Region == "US") %>%
select(-c(Province.State, Country.Region, Long, Lat)) %>%
group_by(status) %>%
summarise_all(funs(sum))## Data being read from JHU/CCSE repository## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~## Reading data from https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv## Data retrieved on 2020-04-26 22:49:33 || Range of dates on data: 2020-01-22--2020-04-26 | Nbr of records: 264## --------------------------------------------------------------------------------## Data being read from JHU/CCSE repository## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~## Reading data from https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv## Data retrieved on 2020-04-26 22:49:34 || Range of dates on data: 2020-01-22--2020-04-26 | Nbr of records: 264## --------------------------------------------------------------------------------## Data being read from JHU/CCSE repository## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~## Reading data from https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv## Data retrieved on 2020-04-26 22:49:34 || Range of dates on data: 2020-01-22--2020-04-26 | Nbr of records: 250## --------------------------------------------------------------------------------ts.us.all = as_tibble(t(as.matrix(ts.us.all)))## Warning: The `x` argument of `as_tibble.matrix()` must have column names if `.name_repair` is omitted as of tibble 2.0.0.
## Using compatibility `.name_repair`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.colnames(ts.us.all) = unlist(ts.us.all[1,])
ts.us.all = ts.us.all[-1,]
ts.us.all$confirmed = as.integer(ts.us.all$confirmed)
ts.us.all$death = as.integer(ts.us.all$death)
ts.us.all$recovered = as.integer(ts.us.all$recovered)
ts.us.all = ts.us.all %>% mutate(active = confirmed - death - recovered)
head(ts.us.all)## # A tibble: 6 x 4
## confirmed death recovered active
## <int> <int> <int> <int>
## 1 1 0 0 1
## 2 1 0 0 1
## 3 2 0 0 2
## 4 2 0 0 2
## 5 5 0 0 5
## 6 5 0 0 52. Exploratory Analysis
2.1 World Map of Confirmed Cases and Deaths
world_map <- map_data("world")
agg.gloabl.map = agg.gloabl
agg.gloabl.map$Country_Region = as.character(agg.gloabl$Country_Region)
agg.gloabl.map$Country_Region[agg.gloabl.map$Country_Region == "US"] = "USA"
agg.gloabl.map$Country_Region[agg.gloabl.map$Country_Region == "United Kingdom"] = "UK"
agg.gloabl.map$Country_Region[agg.gloabl.map$Country_Region == "Korea, South"] = "South Korea"
agg.gloabl.map$Country_Region[agg.gloabl.map$Country_Region == "Taiwan*"] = "Taiwan"
agg.gloabl.map <- inner_join(agg.gloabl.map, world_map, by = c("Country_Region" = "region"))ggplotly(
ggplot(agg.gloabl.map, aes(long, lat, group = group, fill = Confirmed, text = paste0(Country_Region, "\n", "Confirmed Cases:", Confirmed))) +
scale_fill_gradient(low = "gray99", high = "coral3") +
geom_polygon(), tooltip = c("text")
)ggplotly(
ggplot(agg.gloabl.map, aes(long, lat, group = group, fill = log(Confirmed), text = paste0(Country_Region, "\n", "Confirmed Cases:", Confirmed))) +
scale_fill_gradient(low = "gray99", high = "coral3") +
geom_polygon(), tooltip = c("text")
)ggplotly(
ggplot(agg.gloabl.map, aes(long, lat, group = group, fill = Deaths, text = paste0(Country_Region, "\n", "Deaths:", Deaths))) +
scale_fill_gradient(low = "gray99", high = "coral3") +
geom_polygon(), tooltip = c("text")
)ggplotly(
ggplot(agg.gloabl.map, aes(long, lat, group = group, fill = log(Deaths), text = paste0(Country_Region, "\n", "Deaths:", Deaths))) +
scale_fill_gradient(low = "gray99", high = "coral3") +
geom_polygon(), tooltip = c("text")
)2.2 US Map of Confirmed Cases and Deaths
states_map = map_data("state")
agg.us.map = agg.us %>% mutate(Province_State = tolower(Province_State))
agg.us.map = inner_join(agg.us.map, states_map, by = c("Province_State" = "region"))ggplotly(
ggplot(agg.us.map, aes(long, lat, group = group, fill = Confirmed, text = paste0(Province_State, "\n", "Confirmed Cases:", Confirmed))) +
scale_fill_gradient(low = "gray99", high = "coral3") +
geom_polygon(), tooltip = c("text")
)ggplotly(
ggplot(agg.us.map, aes(long, lat, group = group, fill = log(Confirmed), text = paste0(Province_State, "\n", "Confirmed Cases:", Confirmed))) +
scale_fill_gradient(low = "gray99", high = "coral3") +
geom_polygon(), tooltip = c("text")
)ggplotly(
ggplot(agg.us.map, aes(long, lat, group = group, fill = Deaths, text = paste0(sapply(Province_State, Caps), "\n", "Deaths:", Deaths))) +
scale_fill_gradient(low = "gray99", high = "coral3") +
geom_polygon(), tooltip = c("text")
)ggplotly(
ggplot(agg.us.map, aes(long, lat, group = group, fill = log(Deaths), text = paste0(sapply(Province_State, Caps), "\n", "Deaths:", Deaths))) +
scale_fill_gradient(low = "gray99", high = "coral3") +
geom_polygon(), tooltip = c("text")
)2.3 Global Time-Series Plot (Top 10 Countries)
top10_countries = ts.confirmed$Country.Region[1:10]
ts.confirmed.top10 = ts.confirmed %>% filter(Country.Region %in% top10_countries)
df.global.ts = tibble()
n_days = 100
for(i in c(1:10)) {
# get the variables
Count = unlist(ts.confirmed.top10[i,-1])[unlist(ts.confirmed.top10[i,-1]) > n_days]
Day = c(1:length(Count))
Country = rep(top10_countries[i], length(Count))
# append to the pre-specified table
df.new = tibble(Country, Day, Count)
df.global.ts = rbind(df.global.ts, df.new)
}ggplotly(
ggplot(df.global.ts, aes(Day, Count, group = Country, color = Country,
text = paste0(Country, "\n", "Total Cases:", Count))) +
geom_line() +
ggtitle("Global Total Confirmed Cases (Top 10 Countries)") +
xlab("Days After 100 Cases") +
ylab("Total Confirmed Cases"), tooltip = c("text")
)2.4 U.S. Time-Series Plot (Top 10 States)
top10_states = ts.us %>% group_by(State) %>% summarise(max = max(Confirmed)) %>%
arrange(desc(max)) %>% select(State)
top10_states = top10_states$State[1:10]ggplotly(
ggplot(ts.us %>% filter(State %in% top10_states), aes(Date, Confirmed, group = State, color = State,
text = paste0(State, "\n", "Total Cases:", Confirmed))) +
geom_line() +
ggtitle("U.S. Total Confirmed Cases (Top 10 States)") +
xlab("Date") +
ylab("Total Confirmed Cases"), tooltip = c("text")
)2.5 World Trend
ts.world = colSums(ts.confirmed[,-1])
ts.world = tibble("Date" = names(ts.world), "Count" = ts.world)
ts.world$Date = as.Date(ts.world$Date)
ts.world$New = c(NA, diff(ts.world$Count))
coeff = tail(ts.world$Count, 1)/tail(ts.world$New, 1)
newColor = "black"
countColor = "darkorange"
moving_n = 10
ggplot(ts.world, aes(x=Date)) +
geom_line(aes(y=New), color=newColor, alpha = 0.2) +
geom_ma(aes(y=New), ma_fun = SMA, n = moving_n, color=newColor) +
geom_bar(aes(y=Count/coeff), stat="identity", fill=countColor, color=countColor, alpha=.2) +
scale_y_continuous(
# Features of the first axis
name = "Number of New Cases",
# Add a second axis and specify its features
sec.axis = sec_axis(~.*coeff, name="Total Confirmed Cases")
) +
theme(
axis.title.y = element_text(color = newColor, size=11),
axis.title.y.right = element_text(color = countColor, size=11)
)## Warning: Removed 1 row(s) containing missing values (geom_path).
2.6 U.S. Trend
ts.us.confirmed = unlist(ts.confirmed[1, -1])
ts.us.confirmed = tibble("Date" = names(ts.us.confirmed), "Count" = ts.us.confirmed)
ts.us.confirmed$Date = as.Date(ts.us.confirmed$Date)
ts.us.confirmed$New = c(NA, diff(ts.us.confirmed$Count))
coeff = tail(ts.us.confirmed$Count, 1)/tail(ts.us.confirmed$New, 1)
ggplot(ts.us.confirmed, aes(x=Date)) +
geom_line(aes(y=New), color=newColor, alpha = 0.2) +
geom_ma(aes(y=New), ma_fun = SMA, n = moving_n, color=newColor) +
geom_bar(aes(y=Count/coeff), stat="identity", fill=countColor, color=countColor, alpha=.2) +
scale_y_continuous(
# Features of the first axis
name = "Number of New Cases",
# Add a second axis and specify its features
sec.axis = sec_axis(~.*coeff, name="Total Confirmed Cases")
) +
theme(
axis.title.y = element_text(color = newColor, size=11),
axis.title.y.right = element_text(color = countColor, size=11)
)## Warning: Removed 1 row(s) containing missing values (geom_path).
3. U.S. COVID-19 Analysis and Prediction
3.1 U.S. Growth Rate
data = covid19.data("ts-confirmed")## Data being read from JHU/CCSE repository## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~## Reading data from https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv## Data retrieved on 2020-04-26 22:49:53 || Range of dates on data: 2020-01-22--2020-04-26 | Nbr of records: 264## --------------------------------------------------------------------------------growth.rate(data, geo.loc="US")## Processing... US## Loading required package: pheatmap## Loading required package: gplots##
## Attaching package: 'gplots'## The following object is masked from 'package:PerformanceAnalytics':
##
## textplot## The following object is masked from 'package:stats':
##
## lowess
## $Changes
## geo.loc 2020-01-23 2020-01-24 2020-01-25 2020-01-26 2020-01-27 2020-01-28
## 1 US 0 1 0 3 0 0
## 2020-01-29 2020-01-30 2020-01-31 2020-02-01 2020-02-02 2020-02-03 2020-02-04
## 1 0 0 2 1 0 3 0
## 2020-02-05 2020-02-06 2020-02-07 2020-02-08 2020-02-09 2020-02-10 2020-02-11
## 1 0 0 0 0 0 0 1
## 2020-02-12 2020-02-13 2020-02-14 2020-02-15 2020-02-16 2020-02-17 2020-02-18
## 1 0 1 0 0 0 0 0
## 2020-02-19 2020-02-20 2020-02-21 2020-02-22 2020-02-23 2020-02-24 2020-02-25
## 1 0 0 2 0 0 36 0
## 2020-02-26 2020-02-27 2020-02-28 2020-02-29 2020-03-01 2020-03-02 2020-03-03
## 1 6 1 2 8 6 24 20
## 2020-03-04 2020-03-05 2020-03-06 2020-03-07 2020-03-08 2020-03-09 2020-03-10
## 1 31 68 45 140 116 65 376
## 2020-03-11 2020-03-12 2020-03-13 2020-03-14 2020-03-15 2020-03-16 2020-03-17
## 1 322 382 516 548 772 1133 1789
## 2020-03-18 2020-03-19 2020-03-20 2020-03-21 2020-03-22 2020-03-23 2020-03-24
## 1 1362 5964 5526 6327 7676 10567 9893
## 2020-03-25 2020-03-26 2020-03-27 2020-03-28 2020-03-29 2020-03-30 2020-03-31
## 1 12042 18058 17821 19808 19444 20922 26341
## 2020-04-01 2020-04-02 2020-04-03 2020-04-04 2020-04-05 2020-04-06 2020-04-07
## 1 25070 30380 31745 33283 28152 29515 30804
## 2020-04-08 2020-04-09 2020-04-10 2020-04-11 2020-04-12 2020-04-13 2020-04-14
## 1 31533 34126 33755 29861 28917 25306 27051
## 2020-04-15 2020-04-16 2020-04-17 2020-04-18 2020-04-19 2020-04-20 2020-04-21
## 1 28680 31242 32114 32491 26612 25517 27539
## 2020-04-22 2020-04-23 2020-04-24 2020-04-25 2020-04-26
## 1 28486 28819 36188 32796 27631
##
## $Growth.Rate
## geo.loc 2020-01-24 2020-01-25 2020-01-26 2020-01-27 2020-01-28 2020-01-29
## 1 US NA 0 NA 0 NaN NaN
## 2020-01-30 2020-01-31 2020-02-01 2020-02-02 2020-02-03 2020-02-04 2020-02-05
## 1 NaN NA 0.5 0 NA 0 NaN
## 2020-02-06 2020-02-07 2020-02-08 2020-02-09 2020-02-10 2020-02-11 2020-02-12
## 1 NaN NaN NaN NaN NaN NA 0
## 2020-02-13 2020-02-14 2020-02-15 2020-02-16 2020-02-17 2020-02-18 2020-02-19
## 1 NA 0 NaN NaN NaN NaN NaN
## 2020-02-20 2020-02-21 2020-02-22 2020-02-23 2020-02-24 2020-02-25 2020-02-26
## 1 NaN NA 0 NaN NA 0 NA
## 2020-02-27 2020-02-28 2020-02-29 2020-03-01 2020-03-02 2020-03-03 2020-03-04
## 1 0.1666667 2 4 0.75 4 0.8333333 1.55
## 2020-03-05 2020-03-06 2020-03-07 2020-03-08 2020-03-09 2020-03-10 2020-03-11
## 1 2.193548 0.6617647 3.111111 0.8285714 0.5603448 5.784615 0.856383
## 2020-03-12 2020-03-13 2020-03-14 2020-03-15 2020-03-16 2020-03-17 2020-03-18
## 1 1.186335 1.350785 1.062016 1.408759 1.467617 1.578994 0.7613192
## 2020-03-19 2020-03-20 2020-03-21 2020-03-22 2020-03-23 2020-03-24 2020-03-25
## 1 4.378855 0.9265594 1.144951 1.213213 1.376628 0.9362165 1.217224
## 2020-03-26 2020-03-27 2020-03-28 2020-03-29 2020-03-30 2020-03-31 2020-04-01
## 1 1.499585 0.9868756 1.111498 0.9816236 1.076013 1.25901 0.9517482
## 2020-04-02 2020-04-03 2020-04-04 2020-04-05 2020-04-06 2020-04-07 2020-04-08
## 1 1.211807 1.044931 1.048449 0.8458372 1.048416 1.043673 1.023666
## 2020-04-09 2020-04-10 2020-04-11 2020-04-12 2020-04-13 2020-04-14 2020-04-15
## 1 1.082231 0.9891285 0.8846393 0.9683869 0.8751254 1.068956 1.06022
## 2020-04-16 2020-04-17 2020-04-18 2020-04-19 2020-04-20 2020-04-21 2020-04-22
## 1 1.089331 1.027911 1.011739 0.8190576 0.9588531 1.079241 1.034388
## 2020-04-23 2020-04-24 2020-04-25 2020-04-26 NA
## 1 1.01169 1.255699 0.9062673 0.8425113 NA3.2 SEIRD Model Simulation in the U.S.
3.2.1 Model Tuning
SEIRD = function(S0, E0, I0, R0, D0, alpha1, alpha2, c1, c2, beta1, gamma1, eta1, n1) {
Out1 = matrix(0, ncol = 5, nrow = n1) # Output container
for(i in 1:n1) {
S0n = S0
I0n = I0
R0n = R0
E0n = E0
D0n = D0
ind1 = ifelse(i > c1, 1, 0) # intervetion
ind2 = ifelse(i > c2, 1, 0) # intervetion
S0 = max(0, S0n - (alpha1 + alpha2*ind1 + log(i)*0.068*alpha2*ind2)*I0n*S0n)
E0 = max(0, E0 + (alpha1 + alpha2*ind1 + log(i)*0.068*alpha2*ind2)*I0n*S0n - gamma1*E0n)
I0 = max(0, I0n + gamma1*E0n - beta1*I0n - eta1*I0n)
R0 = max(0, R0n + beta1*I0n)
D0 = max(0, D0n + eta1*I0n)
Out1[i,1] = S0
Out1[i,2] = I0
Out1[i,3] = R0
Out1[i,4] = E0
Out1[i,5] = D0
}
Out2 = tibble(S = Out1[,1],
I = Out1[,2],
R = Out1[,3],
E = Out1[,4],
D = Out1[,5])
return(Out2)
}
# SEIRD Model Parameter Tuning
S0 = 327200000 # Initial susceptible
I0 = 1 # Initial infected
R0 = 0 # Initial recovered
E0 = 10 # Initial exposed
D0 = 0 # Initial death
alpha1 = 1/860000000 # S - > E
alpha2 = -1/1250000000 # Intervention
beta1 = 1/140 # I -> R
gamma1 = 1/7 # E -> I
eta1 = 1/230 # I -> D
n1 = nrow(ts.us.all)
Out1 = SEIRD(S0, E0, I0, R0, D0, alpha1, alpha2, 70, 82, beta1, gamma1, eta1, n1)
plot(1:n1, Out1$I, col = "orange", main = "Infected", xlab="Days", ylab="Cases")
lines(1:n1, ts.us.all$active, col = "orange")
plot(1:n1, Out1$D, col = "red", main = "Deaths", xlab="Days", ylab="Cases")
lines(1:n1, ts.us.all$death, col = "red")
plot(1:n1, Out1$R, col = "green", main = "Recovered", xlab="Days", ylab="Cases")
lines(1:n1, ts.us.all$recovered, col = "green")
3.2.2 1000 Days Simulation
SEIRD = function(S0, E0, I0, R0, D0, alpha1, alpha2, c1, c2, beta1, gamma1, eta1, n1) {
Out1 = matrix(0, ncol = 5, nrow = n1) # Output container
for(i in 1:n1) {
S0n = S0
I0n = I0
R0n = R0
E0n = E0
D0n = D0
ind1 = ifelse(i > c1, 1, 0) # intervetion
ind2 = ifelse(i > c2, 1, 0) # intervetion
S0 = max(0, S0n - (alpha1 + alpha2*ind1 + log(i)*0.081*alpha2*ind2)*I0n*S0n)
E0 = max(0, E0 + (alpha1 + alpha2*ind1 + log(i)*0.081*alpha2*ind2)*I0n*S0n - gamma1*E0n)
I0 = max(0, I0n + gamma1*E0n - beta1*I0n - eta1*I0n)
R0 = max(0, R0n + beta1*I0n)
D0 = max(0, D0n + eta1*I0n)
Out1[i,1] = S0
Out1[i,2] = I0
Out1[i,3] = R0
Out1[i,4] = E0
Out1[i,5] = D0
}
Out2 = tibble(S = Out1[,1],
I = Out1[,2],
R = Out1[,3],
E = Out1[,4],
D = Out1[,5])
return(Out2)
}
# Simulation for 1000 days
n2 = n1+1000
Out1 = SEIRD(S0, E0, I0, R0, D0, alpha1, alpha2, 70, 82, beta1, gamma1, eta1, n2)
plot(1:n2, Out1$I, col = "orange", main = "SEIRD Simulation", xlab="Days", ylab="Cases", ylim = c(0, 1500000), type = "l")
lines(1:n2, Out1$D, col = "red", xlab="Days", ylab="Cases", type = "l")
lines(1:n2, Out1$R, col = "green", xlab="Days", ylab="Cases", type = "l")